EEG_reduced

library(tidyverse)

## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(patchwork)
library(rmatio)
library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

load("data/behav.RData")
load("data/load_effects_DFR.RData")
load("data/split_groups_info.RData")

source("helper_fxns/load_EEG_data.R")
source("helper_fxns/corr_spectrum_to_indiv_diff.R")
source("helper_fxns/mutate_for_heatmap.R")
source("helper_fxns/split_ERPs_into_groups.R")
source("helper_fxns/create_TC_for_plot.R")
source("helper_fxns/split_spectrum.R")
source("helper_fxns/paired_freq_plot.R")
source("helper_fxns/select_period_average.R")
source("helper_fxns/split_into_groups.R")
source("helper_fxns/prep_split_for_bar_plots.R")
source("helper_fxns/average_electrodes.R")
source("helper_fxns/average_electrodes_spectrum.R")
source("helper_fxns/plot_tertiles.R")
source("helper_fxns/find_inverted_U_spectra.R")

temp <- read.mat("data/EEG/DFR/ERPS_cl15_midOccip_reformatted.mat")
ERPS_times_DFR <- temp[["all_times"]]
temp <- read.mat("data/EEG/DFR/ERSPS_cl15_midOccip_reformatted.mat")
ERSPS_times_DFR <- temp[["all_times"]]
ERSPS_freqs <- temp[["all_freqs"]]


ERPS_midOccip_DFR <- load_EEG_data("DFR", "ERPS_cl15_midOccip")
ERSPS_midOccip_DFR <- load_EEG_data("DFR", "ERSPS_cl15_midOccip")
ERSPS_Oz_DFR <- load_EEG_data("DFR", "ERSPS_Oz")

ERPS_Pz_DFR <- load_EEG_data("DFR", "ERPS_Pz")
ERPS_P1_DFR <- load_EEG_data("DFR", "ERPS_P1")
ERPS_P2_DFR <- load_EEG_data("DFR", "ERPS_P2")
ERPS_CP1_DFR <- load_EEG_data("DFR", "ERPS_CP1")
ERPS_CP2_DFR <- load_EEG_data("DFR", "ERPS_CP2")
ERPS_POz_DFR <- load_EEG_data("DFR", "ERPS_POz")

ERSPS_Fz_DFR <- load_EEG_data("DFR", "ERSPS_Fz")
ERSPS_F1_DFR <- load_EEG_data("DFR", "ERSPS_F1")
ERSPS_F2_DFR <- load_EEG_data("DFR", "ERSPS_F2")
ERSPS_FC1_DFR <- load_EEG_data("DFR", "ERSPS_FC1")
ERSPS_FC2_DFR <- load_EEG_data("DFR", "ERSPS_FC2")
ERSPS_AFz_DFR <- load_EEG_data("DFR", "ERSPS_AFz")


ERPS_O1_DFR <- load_EEG_data("DFR", "ERPS_O1")
ERPS_O2_DFR <- load_EEG_data("DFR", "ERPS_O2")
ERPS_PO7_DFR <- load_EEG_data("DFR", "ERPS_PO7")
ERPS_PO8_DFR <- load_EEG_data("DFR", "ERPS_PO8")


ERPS_Oavg_DFR <- average_electrodes(list(ERPS_O1_DFR, ERPS_O2_DFR))
ERPS_POavg_DFR <- average_electrodes(list(ERPS_PO7_DFR, ERPS_PO8_DFR))

ERPS_Pavg_DFR <- average_electrodes(list(ERPS_Pz_DFR,ERPS_P1_DFR, ERPS_P2_DFR,ERPS_POz_DFR,ERPS_CP1_DFR, ERPS_CP2_DFR))
ERSPS_Favg_DFR <- average_electrodes_spectrum(list(ERSPS_Fz_DFR, ERSPS_F1_DFR, ERSPS_F2_DFR, ERSPS_AFz_DFR, ERSPS_FC1_DFR, ERSPS_FC2_DFR))

#save(list=c("CDA","CDA_fMRI","ERPS_times_DFR", "ERPS_times_LCD", "ERSPS_times_DFR", "ERSPS_times_LCD", "ERSPS_freqs", "ERPS_midOccip_LCD", "ERSPS_midOccip_LCD", "ERSPS_Oz_LCD","ERSPS_O2_LCD","ERSPS_O1_LCD", "ERPS_Pz_LCD","ERPS_Fz_LCD", "ERPS_midOccip_DFR", "ERSPS_midOccip_DFR", "ERSPS_Oz_DFR","ERPS_Pz_DFR"), file="data/newEEG_data.RData")

rects <- data.frame(xstart=c(0,5500),xend=c(2500,7000),col = "gray")

If we remove EEG from the span calculation, 10% of the subjects change their span group classification.

check_span_groups <- constructs_fMRI 

check_span_groups<- check_span_groups[order(check_span_groups$omnibus_span_no_DFR_MRI),]
check_span_groups$without_MRI <- "low"
check_span_groups$without_MRI[57] <- "not_incl"
check_span_groups$without_MRI[58:113] <- "med"
check_span_groups$without_MRI[114] <- "not_incl"
check_span_groups$without_MRI[115:170] <- "high"

check_span_groups<- check_span_groups[order(check_span_groups$omnibus_span_no_DFR),]
check_span_groups$without_EEG <- "low"
check_span_groups$without_EEG[57] <- "not_incl"
check_span_groups$without_EEG[58:113] <- "med"
check_span_groups$without_EEG[114] <- "not_incl"
check_span_groups$without_EEG[115:170] <- "high"

sum(check_span_groups$without_MRI != check_span_groups$without_EEG, na.rm=TRUE)

## [1] 17

check_span_groups <- merge(check_span_groups, p200_demographics)
colnames(check_span_groups)[10] <- "level"

WM_groups_no_EEG <- list(high = check_span_groups %>% filter(level == "high"),
                         med = check_span_groups %>% filter(level == "med"), 
                         low = check_span_groups %>% filter(level == "low"))

WM_groups_no_EEG[["all"]] <- rbind(WM_groups_no_EEG[["low"]],WM_groups_no_EEG[["med"]],WM_groups_no_EEG[["high"]] )

ERP Data

avg_ERPs_for_plot <- list(
  midOccip_DFR = data.frame(high_load = apply(ERPS_midOccip_DFR[["high_load"]][,2:1922], 2, mean),
                            low_load = apply(ERPS_midOccip_DFR[["low_load"]][,2:1922], 2, mean),
                            load_effect = apply(ERPS_midOccip_DFR[["load_effect"]][,2:1922], 2, mean),
                            time =ERPS_times_DFR),
  Pavg_DFR = data.frame(high_load = apply(ERPS_Pavg_DFR[["high_load"]][,2:1922], 2, mean),
                        low_load = apply(ERPS_Pavg_DFR[["low_load"]][,2:1922], 2, mean),
                        load_effect = apply(ERPS_Pavg_DFR[["load_effect"]][,2:1922], 2, mean),
                        time =ERPS_times_DFR),
  Oavg_DFR = data.frame(high_load = apply(ERPS_Oavg_DFR[["high_load"]][,2:1922], 2, mean),
                        low_load = apply(ERPS_Oavg_DFR[["low_load"]][,2:1922], 2, mean),
                        load_effect = apply(ERPS_Oavg_DFR[["load_effect"]][,2:1922], 2, mean),
                        time =ERPS_times_DFR), 
  POavg_DFR = data.frame(high_load = apply(ERPS_POavg_DFR[["high_load"]][,2:1922], 2, mean),
                         low_load = apply(ERPS_POavg_DFR[["low_load"]][,2:1922], 2, mean),
                         load_effect = apply(ERPS_POavg_DFR[["load_effect"]][,2:1922], 2, mean),
                         time =ERPS_times_DFR)
  
)

plot_list <- list()

for (cluster in seq.int(1,length(avg_ERPs_for_plot))){
  plot_list[[names(avg_ERPs_for_plot)[cluster]]][["indiv_loads"]] <- ggplot(data = avg_ERPs_for_plot[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("ERP")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  plot_list[[names(avg_ERPs_for_plot)[cluster]]][["LE"]] <- ggplot(data = avg_ERPs_for_plot[[cluster]])+
    geom_line(aes(x=time,y=load_effect))+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    ylab("ERP")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

avg_P <- (avg_ERPs_for_plot[["Pavg_DFR"]]$high_load + avg_ERPs_for_plot[["Pavg_DFR"]]$low_load)/2
avg_O <- (avg_ERPs_for_plot[["Oavg_DFR"]]$high_load + avg_ERPs_for_plot[["Oavg_DFR"]]$low_load)/2
avg_PO <- (avg_ERPs_for_plot[["POavg_DFR"]]$high_load + avg_ERPs_for_plot[["POavg_DFR"]]$low_load)/2

avg_plots <- data.frame(time = avg_ERPs_for_plot[["Pavg_DFR"]]$time, avg_P, avg_O, avg_PO)

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_P))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average P electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_O))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average O electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_PO))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average PO electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

cue_P3_range <- c(which(avg_plots$avg_P == max(avg_plots$avg_P[229:290]))-6 , which(avg_plots$avg_P == max(avg_plots$avg_P[229:290]))+6)


probe_P3_range <- c(which(avg_plots$avg_P == max(avg_plots$avg_P[1664:1680]))-6 , which(avg_plots$avg_P == max(avg_plots$avg_P[1664:1680]))+6)

cue_n170_range <- c(which(avg_plots$avg_O == min(avg_plots$avg_O[129:198]))-6, which(avg_plots$avg_O == min(avg_plots$avg_O[129:198]))+6)

probe_n170_range <- c(which(avg_plots$avg_O == min(avg_plots$avg_O[1550:1677]))-6,which(avg_plots$avg_O == min(avg_plots$avg_O[1550:1677]))+6)

All Subjects

During both the cue and probe, we see significant load effects in the N170 components in the average of the O1 and O2 electrodes and the PO7 and PO8 averages. We also see a load effect in the P3 in the parietal electrode average. Both of these show a larger amplitude in the low load than high load.

plot_list[["Pavg_DFR"]][["indiv_loads"]] + plot_list[["Pavg_DFR"]][["LE"]] +
  plot_annotation(title="P avg - DFR")

plot_list[["Oavg_DFR"]][["indiv_loads"]] + plot_list[["Oavg_DFR"]][["LE"]] +
  plot_annotation(title="O avg - DFR")

plot_list[["POavg_DFR"]][["indiv_loads"]] + plot_list[["POavg_DFR"]][["LE"]] +
  plot_annotation(title="PO avg - DFR")

cue_average_midOccip_n170 <- select_period_average(ERPS_midOccip_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)
cue_average_O_n170 <- select_period_average(ERPS_Oavg_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)
cue_average_PO_n170 <- select_period_average(ERPS_POavg_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)

cue_average_P3 <- select_period_average(ERPS_Pavg_DFR,avg_plots$time[cue_P3_range[1]],avg_plots$time[cue_P3_range[2]],ERPS_times_DFR)

probe_average_midOccip_n170 <- select_period_average(ERPS_midOccip_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)
probe_average_O_n170 <- select_period_average(ERPS_Oavg_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)
probe_average_PO_n170 <- select_period_average(ERPS_POavg_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)

probe_average_P3 <- select_period_average(ERPS_Pavg_DFR,avg_plots$time[probe_P3_range[1]],avg_plots$time[probe_P3_range[2]],ERPS_times_DFR)

t.test(cue_average_O_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_O_n170$load_effect
## t = -5.2605, df = 189, p-value = 3.871e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8283886 -0.3765553
## sample estimates:
## mean of x 
## -0.602472

t.test(cue_average_PO_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_PO_n170$load_effect
## t = -6.9534, df = 189, p-value = 5.642e-11
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0704404 -0.5973194
## sample estimates:
##  mean of x 
## -0.8338799

t.test(cue_average_P3$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_P3$load_effect
## t = -5.1679, df = 189, p-value = 5.988e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6709431 -0.3002420
## sample estimates:
##  mean of x 
## -0.4855926

t.test(probe_average_O_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_O_n170$load_effect
## t = 0.24254, df = 189, p-value = 0.8086
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3504347  0.4486914
## sample estimates:
##  mean of x 
## 0.04912837

t.test(probe_average_PO_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_PO_n170$load_effect
## t = -0.28819, df = 189, p-value = 0.7735
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8147653  0.6070440
## sample estimates:
##  mean of x 
## -0.1038607

t.test(probe_average_P3$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_P3$load_effect
## t = -4.9474, df = 189, p-value = 1.657e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7978863 -0.3430026
## sample estimates:
##  mean of x 
## -0.5704444

split_ERPS_midOccip_DFR <- split_ERPs_into_groups(ERPS_midOccip_DFR, WM_groups_no_EEG, ERPS_times_DFR) 
split_ERPS_Pavg_DFR <- split_ERPs_into_groups(ERPS_Pavg_DFR, WM_groups_no_EEG, ERPS_times_DFR) 
split_ERPS_Oavg_DFR <- split_ERPs_into_groups(ERPS_Oavg_DFR, WM_groups_no_EEG, ERPS_times_DFR) 
split_ERPS_POavg_DFR <- split_ERPs_into_groups(ERPS_POavg_DFR, WM_groups_no_EEG, ERPS_times_DFR) 


ERPs_for_plot <- list(mid_Occip_DFR = split_ERPS_midOccip_DFR,
                      Pavg_DFR = split_ERPS_Pavg_DFR, 
                      Oavg_DFR = split_ERPS_Oavg_DFR, 
                      POavg_DFR = split_ERPS_POavg_DFR)

split_ERPs_plot <- create_TC_for_plot(ERPs_for_plot)

Split by WM groups

ggplot(data = split_ERPs_plot[["Pavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("P avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["Pavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("P avg DFR")+
  theme_classic()

ggplot(data = split_ERPs_plot[["Oavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("O avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["Oavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("O avg DFR")+
  theme_classic()

ggplot(data = split_ERPs_plot[["POavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("PO avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["POavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("PO avg DFR")+
  theme_classic()

cue_P3_DFR_split <- split_into_groups(cue_average_P3, WM_groups_no_EEG)
cue_midOccip_n170_DFR_split <- split_into_groups(cue_average_midOccip_n170, WM_groups_no_EEG)
cue_Oavg_n170_DFR_split <- split_into_groups(cue_average_O_n170, WM_groups_no_EEG)
cue_POavg_n170_DFR_split <- split_into_groups(cue_average_PO_n170, WM_groups_no_EEG)


probe_P3_DFR_split <- split_into_groups(probe_average_P3, WM_groups_no_EEG)
probe_midOccip_n170_DFR_split <- split_into_groups(probe_average_midOccip_n170, WM_groups_no_EEG)
probe_Oavg_n170_DFR_split <- split_into_groups(probe_average_O_n170, WM_groups_no_EEG)
probe_POavg_n170_DFR_split <- split_into_groups(probe_average_PO_n170, WM_groups_no_EEG)


cue_P3_DFR_split_prepped <- prep_split_for_bar_plots(cue_P3_DFR_split)
cue_midOccip_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_midOccip_n170_DFR_split)
cue_Oavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_Oavg_n170_DFR_split)
cue_POavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_POavg_n170_DFR_split)


probe_P3_DFR_split_prepped <- prep_split_for_bar_plots(probe_P3_DFR_split)
probe_midOccip_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_midOccip_n170_DFR_split)
probe_Oavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_Oavg_n170_DFR_split)
probe_POavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_POavg_n170_DFR_split)

WM_to_merge <- WM_groups_no_EEG[["all"]][,c(1,10)]

cue_midOccip_n170_DFR_anova <- merge(cue_average_midOccip_n170,WM_to_merge, by="PTID")
probe_midOccip_n170_DFR_anova <- merge(probe_average_midOccip_n170,WM_to_merge, by="PTID")

cue_Oavg_n170_DFR_anova <- merge(cue_average_O_n170,WM_to_merge, by="PTID")
probe_Oavg_n170_DFR_anova <- merge(probe_average_O_n170,WM_to_merge, by="PTID")

cue_POavg_n170_DFR_anova <- merge(cue_average_PO_n170,WM_to_merge, by="PTID")
probe_POavg_n170_DFR_anova <- merge(probe_average_PO_n170,WM_to_merge, by="PTID")

cue_P3_DFR_anova <- merge(cue_average_P3,WM_to_merge, by="PTID")
probe_P3_DFR_anova <- merge(probe_average_P3,WM_to_merge, by="PTID")

Load Effect

We don’t see any difference between load effects

cue_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_P3_LE

cue_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_midOccip_n170_LE


cue_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_Oavg_n170_LE

cue_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_POavg_n170_LE

probe_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_P3_LE

probe_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_midOccip_n170_LE

probe_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_Oavg_n170_LE

probe_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_POavg_n170_LE

cue_P3_LE + probe_P3_LE+
  plot_annotation(title="P3 Load Effect during DFR")

cue_Oavg_n170_LE + probe_Oavg_n170_LE+
  plot_annotation(title="Oavg n170 Load Effect during DFR")

cue_POavg_n170_LE + probe_POavg_n170_LE+
  plot_annotation(title="POavg n170 Load Effect during DFR")

cue_Oavg_n170_LE.aov <- aov(load_effect ~ level, data = cue_Oavg_n170_DFR_anova)
print(summary(cue_Oavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    3.1   1.549   0.633  0.532
## Residuals   158  386.4   2.446

probe_Oavg_n170_LE.aov <- aov(load_effect ~ level, data = probe_Oavg_n170_DFR_anova)
print(summary(probe_Oavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    2.3   1.145   0.128   0.88
## Residuals   158 1417.0   8.968

cue_POavg_n170_LE.aov <- aov(load_effect ~ level, data = cue_POavg_n170_DFR_anova)
print(summary(cue_POavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    2.8   1.411   0.518  0.597
## Residuals   158  430.3   2.723

probe_POavg_n170_LE.aov <- aov(load_effect ~ level, data = probe_POavg_n170_DFR_anova)
print(summary(probe_POavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     40   19.95    0.69  0.503
## Residuals   158   4573   28.94

cue_P3_LE.aov <- aov(load_effect ~ level, data = cue_P3_DFR_anova)
print(summary(cue_P3_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.44  0.2185   0.125  0.882
## Residuals   158 275.51  1.7437

probe_P3_LE.aov <- aov(load_effect ~ level, data = probe_P3_DFR_anova)
print(summary(probe_P3_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    0.6  0.3098   0.114  0.893
## Residuals   158  430.6  2.7252

ERSPs

temp_diff <- constructs_fMRI[,c(1,8)]

omnibus_span_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

temp_diff <- data.frame(p200_data[p200_data$PTID %in% constructs_fMRI$PTID, c(1,7)])

DFR_L3_acc_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

temp_diff <- data.frame(PTID=constructs_fMRI$PTID, LE = p200_delay_DFR$DFR_Load3_Load1)

DFR_fMRI_LE_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

Average over all subs

avg_ERSPS_midOccip_DFR <- mutate_for_heatmap(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]],c(1,2),mean))
avg_ERSPS_Oz_DFR <- mutate_for_heatmap(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]],c(1,2),mean))
avg_ERSPS_Favg_DFR <- mutate_for_heatmap(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]],c(1,2),mean))


ggplot(data=avg_ERSPS_midOccip_DFR, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic()

ggplot(data=avg_ERSPS_Oz_DFR, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic()

Split by WMC

These plots give a good overview, but we will inspect each of the frequency bands more in depth later.

split_ERSPS_midOccip_DFR <- split_spectrum(ERSPS_midOccip_DFR, WM_groups_no_EEG)
split_ERSPS_Oz_DFR <- split_spectrum(ERSPS_Oz_DFR, WM_groups_no_EEG)
split_ERSPS_Favg_DFR <- split_spectrum(ERSPS_Favg_DFR, WM_groups_no_EEG)

split_ERSPS_midOccip_DFR_for_plot <- list()
split_ERSPS_Oz_DFR_for_plot <- list()
split_ERSPS_Favg_DFR_for_plot <- list()

split_ERSPS_midOccip_DFR_plots <- list()
split_ERSPS_Oz_DFR_plots <- list()
split_ERSPS_Favg_DFR_plots <- list()


for (group in c("high","med","low")){
  split_ERSPS_midOccip_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_midOccip_DFR[[group]][["load_effect"]])
  split_ERSPS_Oz_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_Oz_DFR[[group]][["load_effect"]])
  split_ERSPS_Favg_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_Favg_DFR[[group]][["load_effect"]])
  
  
  
  
  split_ERSPS_midOccip_DFR_plots[[group]] <- ggplot(data=split_ERSPS_midOccip_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in midOccip during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic()
  
  split_ERSPS_Oz_DFR_plots[[group]] <- ggplot(data=split_ERSPS_Oz_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in Oz during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic() 
  
  split_ERSPS_Favg_DFR_plots[[group]] <- ggplot(data=split_ERSPS_Favg_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in Favg during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic() 
}

split_ERSPS_midOccip_DFR_plots[["high"]] 

split_ERSPS_midOccip_DFR_plots[["med"]] 

split_ERSPS_midOccip_DFR_plots[["low"]]

split_ERSPS_Oz_DFR_plots[["high"]]

split_ERSPS_Oz_DFR_plots[["med"]]

split_ERSPS_Oz_DFR_plots[["low"]]

Inverted U shape relationships

Taking the linear correlation would obscure any inverted U shaped relationships - so let’s directly look for them by taking the values where abs(med) > abs(low) and abs(med) > abs(high) [taking absolute values to account for when the meaures are negative to baseline, which we see in EEG but don’t necessarily see in fMRI]

inverted_U_midOccip_ERSPS <- find_inverted_U_spectra(split_ERSPS_midOccip_DFR, apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]],c(1,2),mean))
inverted_U_Oz_ERSPS <- find_inverted_U_spectra(split_ERSPS_Oz_DFR, apply(ERSPS_Oz_DFR[["data"]][["load_effect"]],c(1,2),mean))
inverted_U_Favg_ERSPS <- find_inverted_U_spectra(split_ERSPS_Favg_DFR, apply(ERSPS_Favg_DFR[["data"]][["load_effect"]],c(1,2),mean))

Here, we plot the spectra averaged over all subjects, but masked for where there is an inverted U shape relationship.

ggplot(data=inverted_U_midOccip_ERSPS, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Inverted U relationship with WMC in midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic() 

ggplot(data=inverted_U_Oz_ERSPS, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Inverted U relationship with WMC in Oz during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic() 

Correlation to Omnibus Span

Mid occipital cluster shows a positive correlation with span during the encoding and delay period in the beta and low gamma bands, and a negative correlation in the theta through alpha bands during the entire task, but most prominently right at the beginning of the probe period.

Oz shows a strong positive correlation with alpha power during encoding and what looks to be low gamma across the entire task. There also seems to be a negative correlation with beta across the entire task, with particularly strong negative correlations right at the end of the delay period/beginning of probe period.

ggplot(data=omnibus_span_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=omnibus_span_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=omnibus_span_corr[["ERSPS_midOccip_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=omnibus_span_corr[["ERSPS_Oz_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

Correlation to DFR high load accuracy

ggplot(data=DFR_L3_acc_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_L3_acc_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_L3_acc_corr[["ERSPS_midOccip_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_L3_acc_corr[["ERSPS_Oz_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

Correlation to DFR fMRI load effect

There seems to be a strong correlation between fMRI load effect and power in the alpha and beta band across electrodes, particularly during encoding.

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_midOccip_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_Oz_DFR"]][["thresholded_plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

alpha <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean)))),  
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean))))
)

beta <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean)))),
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean))))
)

low_gamma <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean)))),
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean))))
)

alpha_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)
beta_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(beta[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(beta[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(beta[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)


low_gamma_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)

alpha_plot_list <- list()
beta_plot_list <- list()
low_gamma_plot_list <- list()

for (cluster in seq.int(1,length(alpha))){
  alpha_plot_list[[names(alpha_list)[cluster]]][["indiv_loads"]] <- ggplot(data = alpha_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("alpha power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  alpha_plot_list[[names(alpha_list)[cluster]]][["LE"]] <- ggplot(data = alpha_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("alpha power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

for (cluster in seq.int(1,length(beta))){
  beta_plot_list[[names(beta_list)[cluster]]][["indiv_loads"]] <- ggplot(data = beta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("beta power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  beta_plot_list[[names(beta_list)[cluster]]][["LE"]] <- ggplot(data = beta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("beta power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

for (cluster in seq.int(1,length(low_gamma))){
  low_gamma_plot_list[[names(low_gamma_list)[cluster]]][["indiv_loads"]] <- ggplot(data = low_gamma_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("low_gamma power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  low_gamma_plot_list[[names(low_gamma_list)[cluster]]][["LE"]] <- ggplot(data = low_gamma_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("low_gamma power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

Alpha

All Subjects

When we look at the load effects over all subjects, we see significant load effects in the alpha band during the cue and probe period in Oz, mid Occipital cluster and F average electrodes and during the delay period in the Oz electrode. During the cue and probe periods, we see stronger power in the low load, but we see the opposite effect in the delay period.

alpha_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + alpha_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - alpha during DFR")

alpha_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + alpha_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - alpha during DFR")

alpha_cue_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
alpha_cue_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)

alpha_delay_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
alpha_delay_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)

alpha_probe_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
alpha_probe_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(alpha_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_cue_average_midOccip$load_effect
## t = -8.7699, df = 177, p-value = 1.454e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8545668 -0.5406145
## sample estimates:
##  mean of x 
## -0.6975907

t.test(alpha_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_cue_average_Oz$load_effect
## t = -4.4787, df = 189, p-value = 1.298e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4582411 -0.1780102
## sample estimates:
##  mean of x 
## -0.3181257

t.test(alpha_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_delay_average_midOccip$load_effect
## t = -0.1686, df = 177, p-value = 0.8663
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1825836  0.1538409
## sample estimates:
##   mean of x 
## -0.01437135

t.test(alpha_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_delay_average_Oz$load_effect
## t = 2.3001, df = 189, p-value = 0.02253
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.02393801 0.31225787
## sample estimates:
## mean of x 
## 0.1680979

t.test(alpha_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_probe_average_midOccip$load_effect
## t = -3.6404, df = 177, p-value = 0.0003572
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6577849 -0.1953179
## sample estimates:
##  mean of x 
## -0.4265514

t.test(alpha_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_probe_average_Oz$load_effect
## t = -5.941, df = 189, p-value = 1.338e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7887671 -0.3955406
## sample estimates:
##  mean of x 
## -0.5921539

ERSP_behav_corr_data <- constructs_fMRI[,c(1,8)]
ERSP_behav_corr_data$LE <- p200_delay_DFR$DFR_Load3_Load1
ERSP_behav_corr_data[,4] <- p200_data[p200_data$PTID %in% constructs_fMRI$PTID, 7]



alpha_data <- Reduce(function(x,y) merge(x,y, by="PTID"),
                     list(alpha_cue_average_midOccip[,c(1,4)], alpha_cue_average_Oz[,c(1,4)],
                          alpha_delay_average_Oz[,c(1,4)],alpha_probe_average_midOccip[,c(1,4)],
                          alpha_probe_average_Oz[,c(1,4)]))

## Warning in merge.data.frame(x, y, by = "PTID"): column names 'load_effect.x',
## 'load_effect.y' are duplicated in the result

## Warning in merge.data.frame(x, y, by = "PTID"): column names 'load_effect.x',
## 'load_effect.y' are duplicated in the result

colnames(alpha_data) <- c("PTID", "cue_midOccip", "cue_Oz", "delay_Oz", "probe_midOccip", "probe_Oz")

alpha_behav <- merge(ERSP_behav_corr_data,alpha_data, by = "PTID")

for (col in seq.int(5,9)){
  print(ggplot(data = alpha_behav, aes_string(x="omnibus_span_no_DFR", y = colnames(alpha_behav)[col]))+
          geom_point()+
          stat_smooth(method="lm")+
          theme_classic()+
          theme(aspect.ratio=1))
  print(colnames(alpha_behav)[col])
  print(cor.test(alpha_behav$omnibus_span_no_DFR, alpha_behav[,col]))
  
}

## `geom_smooth()` using formula 'y ~ x'

## [1] "cue_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_behav$omnibus_span_no_DFR and alpha_behav[, col]
## t = -1.0087, df = 151, p-value = 0.3147
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.23741074  0.07787438
## sample estimates:
##         cor 
## -0.08181473

## `geom_smooth()` using formula 'y ~ x'

## [1] "cue_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_behav$omnibus_span_no_DFR and alpha_behav[, col]
## t = 0.8596, df = 151, p-value = 0.3914
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08989109  0.22595850
## sample estimates:
##       cor 
## 0.0697824

## `geom_smooth()` using formula 'y ~ x'

## [1] "delay_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_behav$omnibus_span_no_DFR and alpha_behav[, col]
## t = -1.1521, df = 151, p-value = 0.2511
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.24834326  0.06631691
## sample estimates:
##         cor 
## -0.09334348

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_behav$omnibus_span_no_DFR and alpha_behav[, col]
## t = -3.2149, df = 151, p-value = 0.001596
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.39588170 -0.09837576
## sample estimates:
##        cor 
## -0.2531029

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_behav$omnibus_span_no_DFR and alpha_behav[, col]
## t = -1.7674, df = 151, p-value = 0.07919
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.29438935  0.01669387
## sample estimates:
##        cor 
## -0.1423614

alphas_for_plot <- list()
betas_for_plot <- list()
low_gammas_for_plot <- list()
load_list <- c("low_load", "high_load", "load_effect")
ERSPS_list <- list( 
  ERSPS_midOccip_DFR = ERSPS_midOccip_DFR, 
  ERSPS_Oz_DFR = ERSPS_Oz_DFR)

for (cluster in seq.int(1,2)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- alpha[[cluster]][[load]]
  }
  if (substr(names(alpha)[cluster],nchar(names(alpha)[cluster])-4,nchar(names(alpha)[cluster])-2) == "LCD"){
    alphas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_LCD)
  } else{
    alphas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_DFR)
  }
}

for (cluster in seq.int(1,2)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- beta[[cluster]][[load]]
  }
  if (substr(names(beta)[cluster],nchar(names(beta)[cluster])-4,nchar(names(beta)[cluster])-2) == "LCD"){
    betas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_LCD)
  } else{
    betas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_DFR)
  }
}

for (cluster in seq.int(1,2)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- low_gamma[[cluster]][[load]]
  }
  if (substr(names(low_gamma)[cluster],nchar(names(low_gamma)[cluster])-4,nchar(names(low_gamma)[cluster])-2) == "LCD"){
    low_gammas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_LCD)
  } else{
    low_gammas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups_no_EEG, ERSPS_times_DFR)
  }
}

By WM groups

When we split into the working memory capacity groups, we see differences across capacity in the mid occipital cluster during probe period, where we see a more negative load effect in the high capacity group vs the medium capacity group. We also see a significant difference between the low and high capacity subjects in the raw high load power, where the high capacity subjects show more negative power than the low load subjects.

We also can see an inverted U shape during the delay period in the Oz electrode, though it isn’t significant.

split_alphas_plot <- create_TC_for_plot(alphas_for_plot)
split_betas_plot <- create_TC_for_plot(betas_for_plot)
split_low_gammas_plot <- create_TC_for_plot(low_gammas_for_plot)

midOccip_DFR_alpha <- paired_freq_plot(split_alphas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_alpha <- paired_freq_plot(split_alphas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))

midOccip_DFR_alpha[[1]] + midOccip_DFR_alpha[[2]]+
  plot_annotation(title="midOccip - alpha during DFR")+
  plot_layout(guides="collect")

Oz_DFR_alpha[[1]] + Oz_DFR_alpha[[2]]+
  plot_annotation(title="Oz - alpha during DFR")+
  plot_layout(guides="collect")

alpha_cue_Oz_DFR_split <- split_into_groups(alpha_cue_average_Oz, WM_groups_no_EEG)
alpha_cue_midOccip_DFR_split <- split_into_groups(alpha_cue_average_midOccip, WM_groups_no_EEG)

alpha_delay_Oz_DFR_split <- split_into_groups(alpha_delay_average_Oz, WM_groups_no_EEG)
alpha_delay_midOccip_DFR_split <- split_into_groups(alpha_delay_average_midOccip, WM_groups_no_EEG)

alpha_probe_Oz_DFR_split <- split_into_groups(alpha_probe_average_Oz, WM_groups_no_EEG)
alpha_probe_midOccip_DFR_split <- split_into_groups(alpha_probe_average_midOccip, WM_groups_no_EEG)

alpha_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_cue_Oz_DFR_split)
alpha_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_cue_midOccip_DFR_split)


alpha_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_delay_Oz_DFR_split)
alpha_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_delay_midOccip_DFR_split)


alpha_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_probe_Oz_DFR_split)
alpha_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_probe_midOccip_DFR_split)

alpha_cue_midOccip_DFR_anova <- merge(alpha_cue_average_midOccip,WM_to_merge, by="PTID")
alpha_delay_midOccip_DFR_anova <- merge(alpha_delay_average_midOccip,WM_to_merge, by="PTID")
alpha_probe_midOccip_DFR_anova <- merge(alpha_probe_average_midOccip,WM_to_merge, by="PTID")

alpha_cue_Oz_DFR_anova <- merge(alpha_cue_average_Oz,WM_to_merge, by="PTID")
alpha_delay_Oz_DFR_anova <- merge(alpha_delay_average_Oz,WM_to_merge, by="PTID")
alpha_probe_Oz_DFR_anova <- merge(alpha_probe_average_Oz,WM_to_merge, by="PTID")

alpha_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_Oz_LE

alpha_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_midOccip_LE

alpha_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_Oz_LE

alpha_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_midOccip_LE

alpha_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_Oz_LE

alpha_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_midOccip_LE

alpha_cue_Oz_LE + alpha_delay_Oz_LE + alpha_probe_Oz_LE+
  plot_annotation(title="Oz electrode alpha frequency Load Effect during DFR")

alpha_cue_midOccip_LE + alpha_delay_midOccip_LE + alpha_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster alpha frequency Load Effect during DFR")

alpha_cue_midOccip.aov <- aov(load_effect ~ level, data = alpha_cue_midOccip_DFR_anova)
print("alpha - cue midOccip")

## [1] "alpha - cue midOccip"

print(summary(alpha_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.84  0.9212   0.789  0.456
## Residuals   148 172.79  1.1675

alpha_delay_midOccip.aov <- aov(load_effect ~ level, data = alpha_delay_midOccip_DFR_anova)
print("alpha - delay midOccip")

## [1] "alpha - delay midOccip"

print(summary(alpha_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.68  0.8395   0.633  0.532
## Residuals   148 196.27  1.3262

alpha_probe_midOccip.aov <- aov(load_effect ~ level, data = alpha_probe_midOccip_DFR_anova)
print("alpha - probe midOccip")

## [1] "alpha - probe midOccip"

alpha_cue_Oz.aov <- aov(load_effect ~ level, data = alpha_cue_Oz_DFR_anova)
print("alpha - cue Oz")

## [1] "alpha - cue Oz"

print(summary(alpha_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.08  0.0410   0.044  0.957
## Residuals   158 147.56  0.9339

alpha_delay_Oz.aov <- aov(load_effect ~ level, data = alpha_delay_Oz_DFR_anova)
print("alpha - delay Oz")

## [1] "alpha - delay Oz"

print(summary(alpha_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.13  0.5652   0.583   0.56
## Residuals   158 153.23  0.9698

Beta

All Subjects

Overall, we see differences in load effects during cue and probe in all electrodes (low > high) and in delay during mid occipital and Oz (high > low).

beta_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + beta_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - beta during DFR")

beta_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + beta_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - beta during DFR")

beta_cue_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
beta_cue_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)

beta_delay_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
beta_delay_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)

beta_probe_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
beta_probe_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(beta_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_cue_average_midOccip$load_effect
## t = -6.2428, df = 177, p-value = 3.092e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5050572 -0.2624377
## sample estimates:
##  mean of x 
## -0.3837475

t.test(beta_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_cue_average_Oz$load_effect
## t = -3.7009, df = 189, p-value = 0.0002817
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3785452 -0.1153149
## sample estimates:
## mean of x 
##  -0.24693

t.test(beta_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_delay_average_midOccip$load_effect
## t = 2.1457, df = 177, p-value = 0.03326
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01171767 0.28023834
## sample estimates:
## mean of x 
##  0.145978

t.test(beta_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_delay_average_Oz$load_effect
## t = 6.7658, df = 189, p-value = 1.615e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3191722 0.5818761
## sample estimates:
## mean of x 
## 0.4505242

t.test(beta_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_probe_average_midOccip$load_effect
## t = -5.0094, df = 177, p-value = 1.315e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7269556 -0.3160592
## sample estimates:
##  mean of x 
## -0.5215074

t.test(beta_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_probe_average_Oz$load_effect
## t = -6.1079, df = 189, p-value = 5.63e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7840744 -0.4012606
## sample estimates:
##  mean of x 
## -0.5926675

beta_data <- Reduce(function(x,y) merge(x,y, by="PTID"),
                     list(beta_cue_average_midOccip[,c(1,4)], beta_cue_average_Oz[,c(1,4)],
                          beta_delay_average_midOccip[,c(1,4)],beta_delay_average_Oz[,c(1,4)],
                          beta_probe_average_midOccip[,c(1,4)],
                          beta_probe_average_Oz[,c(1,4)]))

## Warning in merge.data.frame(x, y, by = "PTID"): column names 'load_effect.x',
## 'load_effect.y' are duplicated in the result

## Warning in merge.data.frame(x, y, by = "PTID"): column names 'load_effect.x',
## 'load_effect.y' are duplicated in the result

## Warning in merge.data.frame(x, y, by = "PTID"): column names 'load_effect.x',
## 'load_effect.y', 'load_effect.x', 'load_effect.y' are duplicated in the result

colnames(beta_data) <- c("PTID", "cue_midOccip", "cue_Oz", "delay_midOccip","delay_Oz", "probe_midOccip", "probe_Oz")

beta_behav <- merge(ERSP_behav_corr_data,beta_data, by = "PTID")

for (col in seq.int(5,10)){
  print(ggplot(data = beta_behav, aes_string(x="omnibus_span_no_DFR", y = colnames(beta_behav)[col]))+
          geom_point()+
          stat_smooth(method="lm")+
          theme_classic()+
          theme(aspect.ratio=1))
  print(colnames(beta_behav)[col])
  print(cor.test(beta_behav$omnibus_span_no_DFR, beta_behav[,col]))
  
}

## `geom_smooth()` using formula 'y ~ x'

## [1] "cue_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = -0.048243, df = 151, p-value = 0.9616
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1625028  0.1548487
## sample estimates:
##          cor 
## -0.003925906

## `geom_smooth()` using formula 'y ~ x'

## [1] "cue_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = 0.19107, df = 151, p-value = 0.8487
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1434849  0.1737965
## sample estimates:
##        cor 
## 0.01554717

## `geom_smooth()` using formula 'y ~ x'

## [1] "delay_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = -0.53834, df = 151, p-value = 0.5911
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2010493  0.1157143
## sample estimates:
##         cor 
## -0.04376749

## `geom_smooth()` using formula 'y ~ x'

## [1] "delay_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = -1.7392, df = 151, p-value = 0.08403
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.29231708  0.01896054
## sample estimates:
##        cor 
## -0.1401392

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = -2.5399, df = 151, p-value = 0.0121
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.34985360 -0.04518536
## sample estimates:
##        cor 
## -0.2024122

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  beta_behav$omnibus_span_no_DFR and beta_behav[, col]
## t = -1.962, df = 151, p-value = 0.0516
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.308625746  0.001035015
## sample estimates:
##       cor 
## -0.157669

Split by WMC

There are, however, no differences across groups in beta for either the high load trials or load effects.

midOccip_DFR_beta <- paired_freq_plot(split_betas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_beta <- paired_freq_plot(split_betas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))

midOccip_DFR_beta[[1]] + midOccip_DFR_beta[[2]]+
  plot_annotation(title="midOccip - beta during DFR")+
  plot_layout(guides="collect")

Oz_DFR_beta[[1]] + Oz_DFR_beta[[2]]+
  plot_annotation(title="Oz - beta during DFR")+
  plot_layout(guides="collect")

beta_cue_Oz_DFR_split <- split_into_groups(beta_cue_average_Oz, WM_groups)
beta_cue_midOccip_DFR_split <- split_into_groups(beta_cue_average_midOccip, WM_groups)

beta_delay_Oz_DFR_split <- split_into_groups(beta_delay_average_Oz, WM_groups)
beta_delay_midOccip_DFR_split <- split_into_groups(beta_delay_average_midOccip, WM_groups)

beta_probe_Oz_DFR_split <- split_into_groups(beta_probe_average_Oz, WM_groups)
beta_probe_midOccip_DFR_split <- split_into_groups(beta_probe_average_midOccip, WM_groups)

beta_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_cue_Oz_DFR_split)
beta_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_cue_midOccip_DFR_split)

beta_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_delay_Oz_DFR_split)
beta_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_delay_midOccip_DFR_split)

beta_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_probe_Oz_DFR_split)
beta_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_probe_midOccip_DFR_split)

beta_cue_midOccip_DFR_anova <- merge(beta_cue_average_midOccip,WM_to_merge, by="PTID")
beta_delay_midOccip_DFR_anova <- merge(beta_delay_average_midOccip,WM_to_merge, by="PTID")
beta_probe_midOccip_DFR_anova <- merge(beta_probe_average_midOccip,WM_to_merge, by="PTID")

beta_cue_Oz_DFR_anova <- merge(beta_cue_average_Oz,WM_to_merge, by="PTID")
beta_delay_Oz_DFR_anova <- merge(beta_delay_average_Oz,WM_to_merge, by="PTID")
beta_probe_Oz_DFR_anova <- merge(beta_probe_average_Oz,WM_to_merge, by="PTID")

beta_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_Oz_LE

beta_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_midOccip_LE

beta_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_Oz_LE

beta_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_midOccip_LE

beta_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_Oz_LE

beta_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_midOccip_LE

beta_cue_Oz_LE + beta_delay_Oz_LE + beta_probe_Oz_LE+
  plot_annotation(title="Oz electrode beta frequency Load Effect during DFR")

beta_cue_midOccip_LE + beta_delay_midOccip_LE + beta_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster beta frequency Load Effect during DFR")

beta_cue_midOccip.aov <- aov(load_effect ~ level, data = beta_cue_midOccip_DFR_anova)
print("beta - cue midOccip")

## [1] "beta - cue midOccip"

print(summary(beta_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.13  0.0657     0.1  0.905
## Residuals   148  97.61  0.6596

beta_delay_midOccip.aov <- aov(load_effect ~ level, data = beta_delay_midOccip_DFR_anova)
print("beta - delay midOccip")

## [1] "beta - delay midOccip"

print(summary(beta_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.44  0.2190   0.248   0.78
## Residuals   148 130.48  0.8816

beta_cue_Oz.aov <- aov(load_effect ~ level, data = beta_cue_Oz_DFR_anova)
print("beta - cue Oz")

## [1] "beta - cue Oz"

print(summary(beta_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.91  0.4537   0.581  0.561
## Residuals   158 123.47  0.7815

beta_delay_Oz.aov <- aov(load_effect ~ level, data = beta_delay_Oz_DFR_anova)
print("beta - delay Oz")

## [1] "beta - delay Oz"

print(summary(beta_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.66  0.8323   1.054  0.351
## Residuals   158 124.81  0.7899

Low Gamma

All Subjects

Significant load effects during delay and probe in the mid Occipital cluster and Oz electrode (low > high) and F avg during probe (low > high).

low_gamma_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + low_gamma_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - low_gamma during DFR")

low_gamma_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + low_gamma_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - low_gamma during DFR")

low_gamma_cue_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
low_gamma_cue_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)

low_gamma_delay_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
low_gamma_delay_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)

low_gamma_probe_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
low_gamma_probe_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(low_gamma_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_cue_average_midOccip$load_effect
## t = -0.14199, df = 177, p-value = 0.8872
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.06904871  0.05977946
## sample estimates:
##    mean of x 
## -0.004634626

t.test(low_gamma_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_cue_average_Oz$load_effect
## t = -1.4566, df = 189, p-value = 0.1469
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.21352909  0.03212947
## sample estimates:
##   mean of x 
## -0.09069981

t.test(low_gamma_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_delay_average_midOccip$load_effect
## t = -2.4175, df = 177, p-value = 0.01664
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.19328508 -0.01954772
## sample estimates:
##  mean of x 
## -0.1064164

t.test(low_gamma_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_delay_average_Oz$load_effect
## t = -2.0996, df = 189, p-value = 0.03709
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.295546088 -0.009216663
## sample estimates:
##  mean of x 
## -0.1523814

t.test(low_gamma_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_probe_average_midOccip$load_effect
## t = -3.8094, df = 177, p-value = 0.0001919
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6239510 -0.1980945
## sample estimates:
##  mean of x 
## -0.4110228

t.test(low_gamma_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_probe_average_Oz$load_effect
## t = -3.7986, df = 189, p-value = 0.000196
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6970924 -0.2205602
## sample estimates:
##  mean of x 
## -0.4588263

low_gamma_data <- Reduce(function(x,y) merge(x,y, by="PTID"),
                     list(low_gamma_delay_average_Oz[,c(1,4)],
                          low_gamma_probe_average_midOccip[,c(1,4)],
                          low_gamma_probe_average_Oz[,c(1,4)]))
colnames(low_gamma_data) <- c("PTID", "delay_Oz", "probe_midOccip", "probe_Oz")

low_gamma_behav <- merge(ERSP_behav_corr_data,low_gamma_data, by = "PTID")

for (col in seq.int(5,7)){
  print(ggplot(data = low_gamma_behav, aes_string(x="omnibus_span_no_DFR", y = colnames(low_gamma_behav)[col]))+
          geom_point()+
          stat_smooth(method="lm")+
          theme_classic()+
          theme(aspect.ratio=1))
  print(colnames(low_gamma_behav)[col])
  print(cor.test(low_gamma_behav$omnibus_span_no_DFR, low_gamma_behav[,col]))
  
}

## `geom_smooth()` using formula 'y ~ x'

## [1] "delay_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  low_gamma_behav$omnibus_span_no_DFR and low_gamma_behav[, col]
## t = -1.4414, df = 151, p-value = 0.1515
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.27018710  0.04296866
## sample estimates:
##        cor 
## -0.1165038

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_midOccip"
## 
##  Pearson's product-moment correlation
## 
## data:  low_gamma_behav$omnibus_span_no_DFR and low_gamma_behav[, col]
## t = -2.0497, df = 151, p-value = 0.04213
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.314981818 -0.006005518
## sample estimates:
##        cor 
## -0.1645269

## `geom_smooth()` using formula 'y ~ x'

## [1] "probe_Oz"
## 
##  Pearson's product-moment correlation
## 
## data:  low_gamma_behav$omnibus_span_no_DFR and low_gamma_behav[, col]
## t = -2.2114, df = 151, p-value = 0.02851
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.32661742 -0.01897397
## sample estimates:
##        cor 
## -0.1771188

Split by WMC

In the low gamma band, the only load effect that shows any difference is Oz during cue, with a trend for an inverted U shape relationship.

In the raw power, however, we do see differences between groups in the cue and delay periods in the mid occipcital cluster and Oz electrodes. In the cue and delay for mid occipital and delay for Oz, we see a significant low > high difference, and in cue period for the mid occipital and cue and delay for Oz, we see a med > high difference.

midOccip_DFR_low_gamma <- paired_freq_plot(split_low_gammas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))
Oz_DFR_low_gamma <- paired_freq_plot(split_low_gammas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))

midOccip_DFR_low_gamma[[1]] + midOccip_DFR_low_gamma[[2]]+
  plot_annotation(title="midOccip - low_gamma during DFR")+
  plot_layout(guides="collect")

Oz_DFR_low_gamma[[1]] + Oz_DFR_low_gamma[[2]]+
  plot_annotation(title="Oz - low_gamma during DFR")+
  plot_layout(guides="collect")

low_gamma_cue_Oz_DFR_split <- split_into_groups(low_gamma_cue_average_Oz, WM_groups)
low_gamma_cue_midOccip_DFR_split <- split_into_groups(low_gamma_cue_average_midOccip, WM_groups)

low_gamma_delay_Oz_DFR_split <- split_into_groups(low_gamma_delay_average_Oz, WM_groups)
low_gamma_delay_midOccip_DFR_split <- split_into_groups(low_gamma_delay_average_midOccip, WM_groups)

low_gamma_probe_Oz_DFR_split <- split_into_groups(low_gamma_probe_average_Oz, WM_groups)
low_gamma_probe_midOccip_DFR_split <- split_into_groups(low_gamma_probe_average_midOccip, WM_groups)

low_gamma_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_cue_Oz_DFR_split)
low_gamma_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_cue_midOccip_DFR_split)

low_gamma_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_delay_Oz_DFR_split)
low_gamma_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_delay_midOccip_DFR_split)

low_gamma_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_probe_Oz_DFR_split)
low_gamma_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_probe_midOccip_DFR_split)

low_gamma_cue_midOccip_DFR_anova <- merge(low_gamma_cue_average_midOccip,WM_to_merge, by="PTID")
low_gamma_delay_midOccip_DFR_anova <- merge(low_gamma_delay_average_midOccip,WM_to_merge, by="PTID")
low_gamma_probe_midOccip_DFR_anova <- merge(low_gamma_probe_average_midOccip,WM_to_merge, by="PTID")

low_gamma_cue_Oz_DFR_anova <- merge(low_gamma_cue_average_Oz,WM_to_merge, by="PTID")
low_gamma_delay_Oz_DFR_anova <- merge(low_gamma_delay_average_Oz,WM_to_merge, by="PTID")
low_gamma_probe_Oz_DFR_anova <- merge(low_gamma_probe_average_Oz,WM_to_merge, by="PTID")

low_gamma_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_Oz_LE

low_gamma_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_midOccip_LE

low_gamma_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_Oz_LE

low_gamma_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_midOccip_LE

low_gamma_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_Oz_LE

low_gamma_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_midOccip_LE

low_gamma_cue_Oz_LE + low_gamma_delay_Oz_LE + low_gamma_probe_Oz_LE+
  plot_annotation(title="Oz electrode low_gamma frequency Load Effect during DFR")

low_gamma_cue_midOccip_LE + low_gamma_delay_midOccip_LE + low_gamma_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster low_gamma frequency Load Effect during DFR")

low_gamma_cue_midOccip.aov <- aov(load_effect ~ level, data = low_gamma_cue_midOccip_DFR_anova)
print("low_gamma - cue midOccip")

## [1] "low_gamma - cue midOccip"

print(summary(low_gamma_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2  0.102 0.05115     0.3  0.741
## Residuals   148 25.218 0.17039

low_gamma_delay_midOccip.aov <- aov(load_effect ~ level, data = low_gamma_delay_midOccip_DFR_anova)
print("low_gamma - delay midOccip")

## [1] "low_gamma - delay midOccip"

print(summary(low_gamma_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.08  0.0383   0.118  0.889
## Residuals   148  48.01  0.3244

low_gamma_cue_Oz.aov <- aov(load_effect ~ level, data = low_gamma_cue_Oz_DFR_anova)
print("low_gamma - cue Oz")

## [1] "low_gamma - cue Oz"

print(summary(low_gamma_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.94  0.4686   0.696    0.5
## Residuals   158 106.36  0.6732

low_gamma_delay_Oz.aov <- aov(load_effect ~ level, data = low_gamma_delay_Oz_DFR_anova)
print("low_gamma - delay Oz")

## [1] "low_gamma - delay Oz"

print(summary(low_gamma_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.12  0.5606    0.68  0.508
## Residuals   158 130.31  0.8247

CDA

CDA <- data.frame(read.csv("data/EEG/LCD/ERPS_CDA.txt"))
CDA_fMRI <- CDA[CDA$subID %in% constructs_fMRI$PTID,]
CDA_fMRI[is.nan(CDA_fMRI[,2]),2:11] <- NA
colnames(CDA_fMRI)[1] <- "PTID"
CDA_fMRI$CDA_3_1 <- CDA_fMRI$LCD_CDA_3dots - CDA_fMRI$LCD_CDA_1dots
CDA_fMRI$CDA_5_1 <- CDA_fMRI$LCD_CDA_5dots - CDA_fMRI$LCD_CDA_1dots
CDA_fMRI$CDA_5_3 <- CDA_fMRI$LCD_CDA_5dots - CDA_fMRI$LCD_CDA_3dots

Raw Data

First, let’s take a look at the raw data, just for sanity’s sake.

CDA_fMRI %>% 
  select(c("LCD_CDA_1dots", "LCD_CDA_3dots", "LCD_CDA_5dots")) %>% 
  melt() %>% 
  ggplot(aes(x=variable,y=value))+
  geom_jitter()+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="red")+
  xlab("Level")+
  ylab("CDA")

## No id variables; using all as measure variables

## Warning: Removed 24 rows containing non-finite values (stat_summary).

## Warning: Removed 24 rows containing missing values (geom_point).

CDA_fMRI %>% 
  select(c("CDA_3_1", "CDA_5_1", "CDA_5_3")) %>% 
  melt() %>% 
  ggplot(aes(x=variable,y=value))+
  geom_jitter()+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="red")+
  xlab("Level differences")+
  ylab("CDA Load Effects")

## No id variables; using all as measure variables

## Warning: Removed 24 rows containing non-finite values (stat_summary).

## Warning: Removed 24 rows containing missing values (geom_point).

Relationship of Raw Data with Span

Other studies have shown a relationship between the CDA and capacity as measured by LCD K. Our data, only showed a correlation between the CDA in the low load condition (1 dot) and the LCD K max, with no correlation between omnibus span or DFR load effect.

data_to_plot <- data.frame(CDA_fMRI, p200_delay_DFR, omnibus_span = constructs_fMRI$omnibus_span_no_DFR)
data_to_plot <- merge(data_to_plot,p200_cog)

(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$DFR_Load3_Load1
## t = 0.1352, df = 160, p-value = 0.8926
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1437443  0.1646122
## sample estimates:
##        cor 
## 0.01068802

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$DFR_Load3_Load1
## t = -1.524, df = 160, p-value = 0.1295
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.26885513  0.03522659
## sample estimates:
##        cor 
## -0.1196186

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$DFR_Load3_Load1
## t = 0.92497, df = 160, p-value = 0.3564
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08218942  0.22460001
## sample estimates:
##       cor 
## 0.0729303

(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$omnibus_span
## t = -0.81334, df = 160, p-value = 0.4172
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.21622386  0.09092732
## sample estimates:
##         cor 
## -0.06416781

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$omnibus_span
## t = -1.4802, df = 160, p-value = 0.1408
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.26566128  0.03866141
## sample estimates:
##        cor 
## -0.1162269

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$omnibus_span
## t = -1.3219, df = 160, p-value = 0.1881
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.25406112  0.05107662
## sample estimates:
##        cor 
## -0.1039374

(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$XLCD_K_MAX
## t = -0.48689, df = 159, p-value = 0.627
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1921126  0.1167875
## sample estimates:
##         cor 
## -0.03858431

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$XLCD_K_MAX
## t = -1.2834, df = 159, p-value = 0.2012
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.25198309  0.05427012
## sample estimates:
##        cor 
## -0.1012547

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$XLCD_K_MAX
## t = -2.0777, df = 159, p-value = 0.03935
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.309470392 -0.008106698
## sample estimates:
##        cor 
## -0.1625777

CDA_split <- split_into_groups(CDA_fMRI, WM_groups)
CDA_split_plot_data <- prep_split_for_bar_plots(CDA_split)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "LCD_CDA_5dots") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA")+
  ggtitle("CDA - 5 dots")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "LCD_CDA_3dots") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA")+
  ggtitle("CDA - 3 dots")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "LCD_CDA_1dots") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA")+
  ggtitle("CDA - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_anova <- merge(CDA_fMRI,WM_to_merge, by="PTID")

summary(aov(LCD_CDA_1dots ~ level,data=CDA_anova))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.64  0.3182   1.053  0.351
## Residuals   157  47.44  0.3022               
## 8 observations deleted due to missingness

summary(aov(LCD_CDA_3dots ~ level,data=CDA_anova))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.45  0.7228   1.684  0.189
## Residuals   157  67.37  0.4291               
## 8 observations deleted due to missingness

Relationship of Load Effects with Span

Vogel and Machizawa (2004) had mean K = 2.8. Our average K max = 2.68 (so slightly lower). The only thing that even shows a trend with DFR fMRI load effects is in the load effect between the 3 dot and 1 dot condition.

(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_5_3))+
   geom_point()+
   stat_smooth(method="lm")+
   theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_5_1))+
                       geom_point()+
                       stat_smooth(method="lm")+
                       theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_3_1))+
                                           geom_point()+
                                           stat_smooth(method="lm")+
                                           theme_classic()+
                                           theme(aspect.ratio=1))

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$DFR_Load3_Load1
## t = 1.3562, df = 160, p-value = 0.177
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04838862  0.25658013
## sample estimates:
##       cor 
## 0.1066022

cor.test(data_to_plot$CDA_5_1, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$DFR_Load3_Load1
## t = -0.56408, df = 160, p-value = 0.5735
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1973896  0.1104038
## sample estimates:
##         cor 
## -0.04455009

cor.test(data_to_plot$CDA_3_1, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$DFR_Load3_Load1
## t = -1.9681, df = 160, p-value = 0.05079
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3008041900  0.0004669451
## sample estimates:
##        cor 
## -0.1537396

(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_5_3))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_5_1))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_3_1))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1)
                        )

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$omnibus_span
## t = 0.59677, df = 160, p-value = 0.5515
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1078529  0.1998694
## sample estimates:
##        cor 
## 0.04712628

cor.test(data_to_plot$CDA_5_1, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$omnibus_span
## t = 0.30731, df = 160, p-value = 0.759
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1303964  0.1778170
## sample estimates:
##        cor 
## 0.02428743

cor.test(data_to_plot$CDA_3_1, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$omnibus_span
## t = -0.31616, df = 160, p-value = 0.7523
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1784947  0.1297082
## sample estimates:
##       cor 
## -0.024987

mean(data_to_plot$XLCD_K_MAX,na.rm=TRUE)

## [1] 2.686799

(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_5_3))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic()+
    theme(aspect.ratio=1))+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_5_1))+
                              geom_point()+
                              stat_smooth(method="loess")+
                              theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_3_1))+
                                                  geom_point()+
                                                  stat_smooth(method="loess")+
                                                  theme_classic()) 

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$XLCD_K_MAX
## t = 0.68698, df = 159, p-value = 0.4931
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1011252  0.2073309
## sample estimates:
##        cor 
## 0.05440059

cor.test(data_to_plot$CDA_5_1, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$XLCD_K_MAX
## t = 1.118, df = 159, p-value = 0.2652
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06727587  0.23971818
## sample estimates:
##        cor 
## 0.08831801

cor.test(data_to_plot$CDA_3_1, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$XLCD_K_MAX
## t = 0.37401, df = 159, p-value = 0.7089
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1256028  0.1834815
## sample estimates:
##        cor 
## 0.02964807

ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=omnibus_span))+
  geom_point()+
  stat_smooth(method="lm")+
  theme_classic()

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 3 rows containing non-finite values (stat_smooth).

## Warning: Removed 3 rows containing missing values (geom_point).

cor.test(data_to_plot$omnibus_span, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$omnibus_span and data_to_plot$XLCD_K_MAX
## t = 8.8404, df = 165, p-value = 1.384e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4541707 0.6618156
## sample estimates:
##       cor 
## 0.5669319

Numerically, it looks like there is an inverted U shaped relationship for the 5 dots - 1 dot and 3 dot - 1 dot, but those relationships are not significant.

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_5_3") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("5 dots - 3 dots")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_5_1") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("5 dots - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_3_1") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("3 dots - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_5_3.aov <- aov(CDA_5_3 ~ level, data = CDA_anova)
print("CDA - 5 dots - 3 dots")

## [1] "CDA - 5 dots - 3 dots"

print(summary(CDA_5_3.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.41  0.7053   1.112  0.331
## Residuals   157  99.58  0.6343               
## 8 observations deleted due to missingness

CDA_5_1.aov <- aov(CDA_5_1 ~ level, data = CDA_anova)
print("CDA - 5 dots - 1 dot")

## [1] "CDA - 5 dots - 1 dot"

print(summary(CDA_5_1.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.53   0.263   0.465  0.629
## Residuals   157  88.86   0.566               
## 8 observations deleted due to missingness

CDA_3_1.aov <- aov(CDA_3_1 ~ level, data = CDA_anova)
print("CDA - 3 dots - 1 dot")

## [1] "CDA - 3 dots - 1 dot"

print(summary(CDA_3_1.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.12  1.0607   1.797  0.169
## Residuals   157  92.68  0.5903               
## 8 observations deleted due to missingness

Correlations between EEG and fMRI

p200_delay_DFR$PTID <- p200_delay_avgs$PTID 

fMRI_measures <- Reduce(function(x,y) merge(x,y,by="PTID"), 
                        list(p200_delay_DFR, p200_FFA, p200_cue_avgs))
  
fMRI_measures$CUE_L1_avg <- rowMeans(cbind(fMRI_measures$R_CUE_L1,fMRI_measures$L_CUE_L1))
fMRI_measures$DELAY_L1_avg <- rowMeans(cbind(fMRI_measures$R_DELAY_L1,fMRI_measures$L_DELAY_L1))
fMRI_measures$PROBE_L1_avg <- rowMeans(cbind(fMRI_measures$R_PROBE_L1,fMRI_measures$L_PROBE_L1))

fMRI_measures$CUE_L3_avg <- rowMeans(cbind(fMRI_measures$R_CUE_L3,fMRI_measures$L_CUE_L3))
fMRI_measures$DELAY_L3_avg <- rowMeans(cbind(fMRI_measures$R_DELAY_L3,fMRI_measures$L_DELAY_L3))
fMRI_measures$PROBE_L3_avg <- rowMeans(cbind(fMRI_measures$R_PROBE_L3,fMRI_measures$L_PROBE_L3))

fMRI_measures$CUE_LE_avg <- rowMeans(cbind(fMRI_measures$R_CUE_LE,fMRI_measures$L_CUE_LE))
fMRI_measures$DELAY_LE_avg <- rowMeans(cbind(fMRI_measures$R_DELAY_LE,fMRI_measures$L_DELAY_LE))
fMRI_measures$PROBE_LE_avg <- rowMeans(cbind(fMRI_measures$R_PROBE_LE,fMRI_measures$L_PROBE_LE))

ERP_measures <- Reduce(function(x,y) merge(x, y, by ="PTID"), 
                       list(cue_average_midOccip_n170, cue_average_P3, 
                            probe_average_midOccip_n170, probe_average_P3))

colnames(ERP_measures) <- c("PTID", "low_load_cue_midOccip_n170", "high_load_cue_midOccip_n170", "load_effect_cue_midOccip_n170", "low_load_cue_P3", "high_load_cue_P3", "load_effect_cue_P3", "low_load_probe_midOccip_n170", "high_load_probe_midOccip_n170", "load_effect_probe_midOccip_n170", "low_load_probe_P3", "high_load_probe_P3", "load_effect_probe_P3")

ERSP_measures_alpha <- Reduce(function(x,y) merge(x,y,by="PTID"), 
                              list(alpha_cue_average_midOccip,
                                   alpha_delay_average_midOccip,
                                   alpha_cue_average_Oz, 
                                   alpha_delay_average_Oz, 
                                   alpha_probe_average_midOccip, 
                                   alpha_probe_average_Oz))


colnames(ERSP_measures_alpha) <- c("PTID","high_alpha_cue_midOccip", "low_alpha_cue_midOccip", "LE_alpha_cue_midOccip", "high_alpha_delay_midOccip", "low_alpha_delay_midOccip", "LE_alpha_delay_midOccip", "high_alpha_cue_Oz", "low_alpha_cue_Oz", "LE_alpha_cue_Oz", "high_alpha_delay_Oz", "low_alpha_delay_Oz", "LE_alpha_delay_Oz", "high_alpha_probe_midOccip", "low_alpha_probe_midOccip", "LE_alpha_probe_midOccip","high_alpha_probe_Oz", "low_alpha_probe_Oz", "LE_alpha_probe_Oz" )

ERSP_measures_beta <- Reduce(function(x,y) merge(x,y,by="PTID"), 
                              list(beta_cue_average_midOccip,
                                   beta_delay_average_midOccip,
                                   beta_cue_average_Oz, 
                                   beta_delay_average_Oz, 
                                   beta_probe_average_midOccip, 
                                   beta_probe_average_Oz))


colnames(ERSP_measures_beta) <- c("PTID","high_beta_cue_midOccip", "low_beta_cue_midOccip", "LE_beta_cue_midOccip", "high_beta_delay_midOccip", "low_beta_delay_midOccip", "LE_beta_delay_midOccip", "high_beta_cue_Oz", "low_beta_cue_Oz", "LE_beta_cue_Oz", "high_beta_delay_Oz", "low_beta_delay_Oz", "LE_beta_delay_Oz", "high_beta_probe_midOccip", "low_beta_probe_midOccip", "LE_beta_probe_midOccip","high_beta_probe_Oz", "low_beta_probe_Oz", "LE_beta_probe_Oz" )

ERSP_measures <- merge(ERSP_measures_alpha, ERSP_measures_beta, by= "PTID")

ERPs

ERP_fMRI_plots <- list()

ERP_fMRI <- merge(fMRI_measures, ERP_measures, by= "PTID")
sig_linear_ERP_fMRI <- list()

cumulative_df_ERP <- data.frame(matrix(nrow=1,ncol=4))
colnames(cumulative_df_ERP) <- c("ERP_measure", "fMRI_measure", "corr", "p_val")

plot_count <- 1
sig_plot_count <- 1
for (ERP in c(36:45)){
  for (fMRI in c(2:4,23,25:33)){
    ERP_fMRI_plots[[plot_count]] <- ggplot(ERP_fMRI, aes_string(y=colnames(ERP_fMRI)[ERP],x=colnames(ERP_fMRI)[fMRI]))+
      geom_point()+
      stat_smooth(method="lm")+
      theme_classic()+
      theme(aspect.ratio=1)
    
    cor_results <- cor.test(ERP_fMRI[,ERP], ERP_fMRI[,fMRI])
    if (cor_results$p.value < 0.05){
      sig_linear_ERP_fMRI[[sig_plot_count]] <- ERP_fMRI_plots[[plot_count]]
      sig_plot_count <- sig_plot_count +1 
      cumulative_df_ERP <- rbind(cumulative_df_ERP, c(colnames(ERP_fMRI)[ERP], colnames(ERP_fMRI)[fMRI], cor_results$estimate, cor_results$p.value))
    } 
    plot_count <- plot_count+1
  }
}

ERP vs fMRI

It’s hard to interpret these since none of them are comparing like, but just wanted to leave it here for the record.

cumulative_df_ERP <- cumulative_df_ERP[2:nrow(cumulative_df_ERP),]

cumulative_df_ERP %>% 
  kable(format = "html", escape = F) %>%
  kable_styling("striped", full_width = F) 

	ERP_measure	fMRI_measure	corr	p_val
2	load_effect_cue_midOccip_n170	CUE_L1_avg	0.168360585458305	0.0374979529251927
3	low_load_cue_P3	CUE_LE_avg	0.182757533804539	0.0237514841131259
4	load_effect_cue_P3	DELAY_L3_avg	0.209448914546025	0.00936647450051395
5	high_load_probe_P3	DFR_Load3_Load1	-0.17219145504146	0.0333082803235829
6	high_load_probe_P3	DELAY_L1_avg	0.180678759270006	0.0254191873879704
7	load_effect_probe_P3	DELAY_L1_avg	0.186714147764717	0.0208357542948427
8	load_effect_probe_P3	CUE_L3_avg	-0.201741313033054	0.0123946284661481
9	load_effect_probe_P3	PROBE_L3_avg	-0.183907436137367	0.0228699996014812
10	load_effect_probe_P3	CUE_LE_avg	-0.203552413827566	0.011614854390257

for (idx in seq.int(1,length(sig_linear_ERP_fMRI))){
  print(sig_linear_ERP_fMRI[[idx]])
}

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

ERSP vs fMRI

ERSP_fMRI_plots <- list()

ERSP_fMRI <- merge(fMRI_measures, ERSP_measures, by= "PTID")
sig_linear_ERSP_fMRI <- list()

cumulative_df_ERSP <- data.frame(matrix(nrow=1,ncol=4))
colnames(cumulative_df_ERSP) <- c("ERSP_measure", "fMRI_measure", "corr", "p_value")

plot_count <- 1
sig_plot_count <- 1
for (ERSP in seq.int(34,69)){
  for (fMRI in c(2:4,23,25:33)){
    ERSP_fMRI_plots[[plot_count]] <- ggplot(ERSP_fMRI, aes_string(y=colnames(ERSP_fMRI)[ERSP],x=colnames(ERSP_fMRI)[fMRI]))+
      geom_point()+
      stat_smooth(method="lm")+
      theme_classic()+
      theme(aspect.ratio=1)
    
    cor_results <- cor.test(ERSP_fMRI[,ERSP], ERSP_fMRI[,fMRI])
    if (cor_results$p.value < 0.05){
      sig_linear_ERSP_fMRI[[sig_plot_count]] <- ERSP_fMRI_plots[[plot_count]]
      sig_plot_count <- sig_plot_count +1 
      cumulative_df_ERSP <- rbind(cumulative_df_ERSP, c(colnames(ERSP_fMRI)[ERSP], colnames(ERSP_fMRI)[fMRI], cor_results$estimate, cor_results$p.value))
      
    } 
    plot_count <- plot_count+1
  }
}

Linear Correlations

Again, hard to interpret because none of the comparisons are like to like, but putting them here.

cumulative_df_ERSP <- cumulative_df_ERSP[2:nrow(cumulative_df_ERSP),]

cumulative_df_ERSP %>% 
  kable(format = "html", escape = F) %>%
  kable_styling("striped", full_width = F) 

	ERSP_measure	fMRI_measure	corr	p_value
2	LE_alpha_delay_midOccip	PROBE_LE_avg	-0.200627928494351	0.0128964603748167
3	high_alpha_cue_Oz	DELAY_LE_avg	-0.18747120252265	0.0203145529047908
4	low_alpha_cue_Oz	DELAY_LE_avg	-0.199602098790874	0.0133744910905956
5	LE_alpha_delay_Oz	PROBE_L1_avg	0.160372075852926	0.0476741304308567
6	LE_alpha_delay_Oz	PROBE_LE_avg	-0.274928099590905	0.000583209202042236
7	high_alpha_probe_Oz	PROBE_L1_avg	0.187015738190036	0.0206267548072726
8	low_alpha_probe_Oz	PROBE_L3_avg	0.162744819578852	0.0444360412475263
9	low_beta_delay_midOccip	CUE_L1_avg	-0.160167982288225	0.0479616016647425
10	low_beta_delay_midOccip	PROBE_L1_avg	-0.161441016959498	0.0461918234128772
11	high_beta_cue_Oz	DELAY_L1_avg	0.217068876209105	0.00703488110260297
12	high_beta_cue_Oz	DELAY_LE_avg	-0.200498263090963	0.0129560455101394
13	low_beta_cue_Oz	DELAY_L1_avg	0.190270081226475	0.0184842091585648
14	low_beta_cue_Oz	DELAY_LE_avg	-0.264271989955132	0.000963603621348035
15	LE_beta_cue_Oz	CUE_L3_avg	-0.165049035878733	0.0414693924103712
16	LE_beta_delay_Oz	PROBE_LE_avg	-0.204706181839383	0.0111408706407856

for (idx in seq.int(1,length(sig_linear_ERSP_fMRI))){
  print(sig_linear_ERSP_fMRI[[idx]])
}

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'